1. ARMA-GARCH model

Let $Y_t$ be the observatoin. (Log-difference of the financial time series) Then, \[ Y_t \sim \mbox{ARMA}(p,q)\\\\ Y_t = \phi_1 Y_{t-1} + e_t + \theta_1 e_{t-1} \] with GARCH error \[ e_t \hspace3mm = \hspace3mm \sigma_t \epsilon_t \hspace{5mm} \epsilon_t \sim N(0,1)\\\\ \sigma_t^2 \hspace3mm = \hspace3mm \omega + \alpha Y_t^2 + \beta \sigma_t^2 \]

$Y_t$ is ARMA with error $\epsilon_t$. $\epsilon_t$ is GARCH. \[ \begin{align} \Phi(B) Y_t &= \Theta(B) \epsilon_t \\ \\ \epsilon_t &= \sigma_t e_t \hspace10mm e_t \sim_{iid} N(0,1) \\\\ \sigma_t^2 &= \omega + \alpha \epsilon_{t-1}^2 + \beta \sigma^2_{t-1} \end{align} \]

Ex: Daily SPY

Daily Price of SP500 ETF (SPY) from Jan 02 2000 to Dec 31 2014

  D <- read.csv("http://gozips.uakron.edu/~nmimoto/pages/datasets/SPY.csv", header=T)

## Warning in if (!header) rlabp <- FALSE: the condition has length > 1 and
## only the first element will be used

## Warning in if (header) {: the condition has length > 1 and only the first
## element will be used

  source("http://gozips.uakron.edu/~nmimoto/477/TS_R-90.txt")
    
  head(D)

##   ixDay      Date   Weekday X SPY.Open SPY.High SPY.Low SPY.Close
## 1     1  1/3/2000    Monday 1   148.25   148.25  143.88    145.44
## 2     2  1/4/2000   Tuesday 2   143.53   144.06  139.64    139.75
## 3     3  1/5/2000 Wednesday 3   139.94   141.53  137.25    140.00
## 4     4  1/6/2000  Thursday 4   139.62   141.50  137.75    137.75
## 5     5  1/7/2000    Friday 5   140.31   145.75  140.06    145.75
## 6     6 1/10/2000    Monday 6   146.25   146.91  145.03    146.25
##   SPY.Volume SPY.Adjusted
## 1    8164300       110.33
## 2    8089800       106.01
## 3   12177900       106.20
## 4    6227200       104.50
## 5    8066500       110.57
## 6    5741700       110.94

  X <- ts(D[,7], start=1)
  plot(X)

  X1 <- diff(log(X))
  plot(X1)

  layout(matrix(1:2, 1, 2))
  acf(X1)
  acf(X1^2)

  layout(1)


  # Fit GARCH model 
  library(fGarch)
  Fit1 <-  garchFit(~ garch(1,1), data=X1, cond.dist="norm",  include.mean = FALSE, trace = FALSE)
    
  ##  Fit1@fit$par          #-- estimated parameters
  ##  Fit1@residuals        #-- this is not the garch residuals!!!! same as X1
  ##  Fit1@sigma.t          #-- estimated sig_t
  ##  Fit1@residuals/Fit1@sigma.t       #-- this is the (standardized) GARCH residuals

  print(Fit1@fit$ics)   #-- AIC and BIC are here

##       AIC       BIC       SIC      HQIC 
## -6.470319 -6.465360 -6.470321 -6.468556

  plot(X1)

  plot(Fit1@sigma.t, type="l")   #- estimated Sigma_t

  plot(X1)
  lines( 1.96*Fit1@sigma.t, type="l", col="red")
  lines(-1.96*Fit1@sigma.t, type="l", col="red")

  res1  <- X1/Fit1@sigma.t   #-- this is the residuals
  Randomness.tests(res1)

##   B-L test H0: the sereis is uncorrelated
##   M-L test H0: the square of the sereis is uncorrelated
##   J-B test H0: the sereis came from Normal distribution
##   SD         : Standard Deviation of the series

##      BL15 BL20 BL25  ML15  ML20 JB    SD
## [1,]    0    0    0 0.719 0.855  0 1.001

2. Example: SPY

  library(xts)

  #--- ARMA-GARCH estimation example ---
  D <- read.csv("http://gozips.uakron.edu/~nmimoto/pages/datasets/SPY.csv", header=T)

## Warning in if (!header) rlabp <- FALSE: the condition has length > 1 and
## only the first element will be used

## Warning in if (header) {: the condition has length > 1 and only the first
## element will be used

  head(D)

##   ixDay      Date   Weekday X SPY.Open SPY.High SPY.Low SPY.Close
## 1     1  1/3/2000    Monday 1   148.25   148.25  143.88    145.44
## 2     2  1/4/2000   Tuesday 2   143.53   144.06  139.64    139.75
## 3     3  1/5/2000 Wednesday 3   139.94   141.53  137.25    140.00
## 4     4  1/6/2000  Thursday 4   139.62   141.50  137.75    137.75
## 5     5  1/7/2000    Friday 5   140.31   145.75  140.06    145.75
## 6     6 1/10/2000    Monday 6   146.25   146.91  145.03    146.25
##   SPY.Volume SPY.Adjusted
## 1    8164300       110.33
## 2    8089800       106.01
## 3   12177900       106.20
## 4    6227200       104.50
## 5    8066500       110.57
## 6    5741700       110.94

  X <- xts(D[,7], order.by=as.Date(D[,"Date"], format="%m/%d/%Y"))    #- X is a xts object
  plot(X)

  X1 <- diff(log(X))[-1]
  plot(X1)

  layout(matrix(1:2,1,2))
  acf(X1)
  acf(X1^2)

  layout(1)

  Randomness.tests(as.numeric(X1))

##   B-L test H0: the sereis is uncorrelated
##   M-L test H0: the square of the sereis is uncorrelated
##   J-B test H0: the sereis came from Normal distribution
##   SD         : Standard Deviation of the series

##      BL15 BL20 BL25 ML15 ML20 JB    SD
## [1,]    0    0    0    0    0  0 0.012

  library(fGarch)
  
 
  #- Estimae parameters of ARMA and GARCH at the same time 
  Fit1 <-  garchFit(~ arma(1,1) + garch(1,1), data=X1, cond.dist="norm",  include.mean = FALSE, trace = FALSE)
  summary(Fit1)

## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~arma(1, 1) + garch(1, 1), data = X1, cond.dist = "norm", 
##     include.mean = FALSE, trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ arma(1, 1) + garch(1, 1)
## <environment: 0x7f8ff5b69390>
##  [data = X1]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##         ar1          ma1        omega       alpha1        beta1  
## -1.4401e-01   2.4903e-01   1.5994e-06   1.0580e-01   8.8436e-01  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## ar1    -1.440e-01   1.317e-01   -1.093   0.2743    
## ma1     2.490e-01   1.302e-01    1.913   0.0557 .  
## omega   1.599e-06   2.683e-07    5.961 2.51e-09 ***
## alpha1  1.058e-01   9.451e-03   11.195  < 2e-16 ***
## beta1   8.844e-01   9.345e-03   94.635  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  12225.18    normalized:  3.241034 
## 
## Description:
##  Wed Apr 11 14:54:21 2018 by user:  
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value  
##  Jarque-Bera Test   R    Chi^2  2177.437  0        
##  Shapiro-Wilk Test  R    W      0.9765013 0        
##  Ljung-Box Test     R    Q(10)  11.66489  0.3081103
##  Ljung-Box Test     R    Q(15)  16.35115  0.3591053
##  Ljung-Box Test     R    Q(20)  23.32912  0.2729244
##  Ljung-Box Test     R^2  Q(10)  8.715104  0.5593378
##  Ljung-Box Test     R^2  Q(15)  11.82243  0.692417 
##  Ljung-Box Test     R^2  Q(20)  13.78888  0.8410369
##  LM Arch Test       R    TR^2   8.749246  0.7241858
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -6.479417 -6.471151 -6.479420 -6.476478

  ##Fit1@fit$par          #-- estimated parameters
  ##Fit1@residuals        #-- this is ARCH residuals. (before GARCH)
  ##print(Fit1@fit$ics)   #-- AIC and BIC are here
  ##Randomness.tests(res1)
  ##plot(Fit1)            #-- choose option 13 for residual qq plot

  Fit2 <-  garchFit(~ arma(0,1) + garch(1,1), data=X1, cond.dist="norm",  include.mean = FALSE, trace = FALSE)
  summary(Fit2)

## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~arma(0, 1) + garch(1, 1), data = X1, cond.dist = "norm", 
##     include.mean = FALSE, trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ arma(0, 1) + garch(1, 1)
## <environment: 0x7f8ff8faae68>
##  [data = X1]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##        ma1       omega      alpha1       beta1  
## 1.0525e-01  1.5932e-06  1.0567e-01  8.8453e-01  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## ma1    1.053e-01   1.808e-02    5.821 5.84e-09 ***
## omega  1.593e-06   2.676e-07    5.953 2.63e-09 ***
## alpha1 1.057e-01   9.442e-03   11.191  < 2e-16 ***
## beta1  8.845e-01   9.333e-03   94.771  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  12224.63    normalized:  3.240888 
## 
## Description:
##  Wed Apr 11 14:54:22 2018 by user:  
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value  
##  Jarque-Bera Test   R    Chi^2  2127.103  0        
##  Shapiro-Wilk Test  R    W      0.9767988 0        
##  Ljung-Box Test     R    Q(10)  14.02298  0.1719459
##  Ljung-Box Test     R    Q(15)  18.63737  0.2306394
##  Ljung-Box Test     R    Q(20)  25.79804  0.1725862
##  Ljung-Box Test     R^2  Q(10)  8.686103  0.56213  
##  Ljung-Box Test     R^2  Q(15)  11.83827  0.6912271
##  Ljung-Box Test     R^2  Q(20)  13.80512  0.8402378
##  LM Arch Test       R    TR^2   8.718929  0.7267317
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -6.479655 -6.473043 -6.479657 -6.477304

  Fit3 <-  garchFit(~ arma(0,2) + garch(1,1), data=X1, cond.dist="norm",  include.mean = FALSE, trace = FALSE)
  summary(Fit3)

## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~arma(0, 2) + garch(1, 1), data = X1, cond.dist = "norm", 
##     include.mean = FALSE, trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ arma(0, 2) + garch(1, 1)
## <environment: 0x7f8ff2f6c548>
##  [data = X1]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##         ma1          ma2        omega       alpha1        beta1  
##  1.0242e-01  -2.7523e-02   1.5971e-06   1.0548e-01   8.8469e-01  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## ma1     1.024e-01   1.771e-02    5.782 7.37e-09 ***
## ma2    -2.752e-02   1.761e-02   -1.563    0.118    
## omega   1.597e-06   2.681e-07    5.958 2.56e-09 ***
## alpha1  1.055e-01   9.429e-03   11.187  < 2e-16 ***
## beta1   8.847e-01   9.326e-03   94.861  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  12226.29    normalized:  3.241329 
## 
## Description:
##  Wed Apr 11 14:54:22 2018 by user:  
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value  
##  Jarque-Bera Test   R    Chi^2  2228.007  0        
##  Shapiro-Wilk Test  R    W      0.9761852 0        
##  Ljung-Box Test     R    Q(10)  10.67642  0.3832771
##  Ljung-Box Test     R    Q(15)  15.33805  0.4273523
##  Ljung-Box Test     R    Q(20)  21.88992  0.3465134
##  Ljung-Box Test     R^2  Q(10)  8.906359  0.541013 
##  Ljung-Box Test     R^2  Q(15)  11.88314  0.6878524
##  Ljung-Box Test     R^2  Q(20)  13.86514  0.8372674
##  LM Arch Test       R    TR^2   8.853235  0.7154123
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -6.480008 -6.471742 -6.480011 -6.477069

  sig.t <- xts(Fit2@sigma.t,                order.by=index(X1))  #-- estimated sig_t fo GARCH
  res1  <- xts(Fit2@residuals/Fit2@sigma.t, order.by=index(X1))   #-- this is the (standardized) ARMA-GARCH residuals

  Randomness.tests(as.numeric(res1))

##   B-L test H0: the sereis is uncorrelated
##   M-L test H0: the square of the sereis is uncorrelated
##   J-B test H0: the sereis came from Normal distribution
##   SD         : Standard Deviation of the series

##       BL15  BL20  BL25  ML15 ML20 JB    SD
## [1,] 0.231 0.173 0.172 0.691 0.84  0 1.001

  #- plot log-return and sigma_t -
  layout(matrix(1:2, 2, 1, byrow=T))
  plot(X1, main="SPY daily log-return")
  plot(sig.t, type="l", main="estimated sigma_t")   #- estimated Sigma_t

  #- plot log-return overlay with sigma_t -
  layout(1)
  plot(X1)

  lines( 1.96*sig.t, type="l", col="red")

  lines(-1.96*sig.t, type="l", col="red")

  #--- h-step ahead prediction from ARMA-GARCH ---
  h = 15    #- predict 15 days ahead

  X.pred  <- xts(predict(Fit2, n.ahead=h)[,1], order.by=index(X1)[length(X1)]+(1:h))
  SD.pred <- xts(predict(Fit2, n.ahead=h)[,3], order.by=index(X1)[length(X1)]+(1:h))

  #- plot log-return overlay with sigma_t -
  plot(rbind(X1["2014"], X.pred))
  lines(X.pred, col="red")
  lines( 1.96*SD.pred, type="l", col="red")
  lines(-1.96*SD.pred, type="l", col="red")

  #--- Rolling 1-step prediction from ARMA-GARCH
  Y <- X1
    
  Y.pred1 <- sig.pred1 <- date.stp <- numeric(0)
  for (i in 1:1000){
    T <- Y[(1:250)+(i-1)]
    Fit44 <-  garchFit(~ arma(0,1) + garch(1,1), data=T, cond.dist="norm", include.mean = FALSE, trace = FALSE)
    Y.pred1[i]   <- predict(Fit44, n.ahead=1)[,1]
    sig.pred1[i] <- predict(Fit44, n.ahead=1)[,3]
    date.stp[i]  <- index(Y[250+i])
  }

## Warning in sqrt(diag(fit$cvar)): NaNs produced

  Y.pred   <- xts(Y.pred1,   order.by=as.Date(date.stp))
  sig.pred <- xts(sig.pred1, order.by=as.Date(date.stp))


  #- Compare predicted vs actual -
  plot(Y["2000::2004"]);   lines(Y.pred, col="red");   lines( 1.96*sig.pred, col="red");   lines(-1.96*sig.pred, col="red")

  #- Actual vs Predicted
  head(  cbind(Y[index(Y.pred)], Y.pred)  )

##                      ..1           ..2
## 2000-12-29 -0.0106150060  0.0009925981
## 2001-01-02 -0.0280598120 -0.0012141588
## 2001-01-03  0.0007836377 -0.0029252764
## 2001-01-04  0.0409786484  0.0003025638
## 2001-01-05 -0.0290649392  0.0040535907
## 2001-01-08 -0.0116787375 -0.0026158966

  plot(Y[index(Y.pred)] - Y.pred)

  mean( sign(Y[index(Y.pred)])==sign(Y.pred) )  #- Can Y.pred guess up/down ?

## [1] 0.53

  mean( sign(Y[index(Y.pred)])==1 )             #- num of days SPY went up

## [1] 0.508

3. Cheking Conditional Distribution

Empirical Distribution Function

Sample version of CDF \[ \hat F_n(x) \hspace3mm = \hspace3mm \frac{1}{n} \sum_{i=1}^n I(x_i < x) \]

  X22 <- rnorm(10)
  X22

##  [1]  0.33979142 -0.28412857 -0.15940948 -0.07124178  0.43093553
##  [6] -0.30531458  1.89158391  1.57043769  1.26749188 -0.38981266

  plot(ecdf(X22))

Kolmogolov-Smirnov Test

In iid setting, K-S test can be used to test hypothesis $X_i \sim F$, and

\[ H_0: F = F_0 \hspace5mm vs \hspace5mm H_A: F \ne F_0 \hspace10mm \mbox{{\scriptsize ($F_0$ is completely specified) } } \]
K-S test uses Empirical Distibution (EDF, ECDF), and calculate the test statistic \[ KS \hspace3mm = \hspace3mm \sup_x |\hat F_n(x) - F(x) | \]
Where is $(\frac{i}{n}, X_{(i)})$ and $(\frac{i-1}{n}, X_{(i)})$?

\[ \begin{align} D_n &= \sup_x |\hat F_n(x) - F_0(x)| \\\\ &= \max_i \Big( \Big|\frac{i}{n} - F_0(x_{(i)}) \Big| , \hspace3mm \Big|\frac{i-1}{n} - F_0(x_{(i)}) \Big| \Big) \\ \\ \end{align} \]
If $X_i$ are indeed from $F_0$, \[ D_n \hspace3mm =_d \hspace3mm \max_i \Big( \Big|\frac{i}{n} - U_{(i)} \Big| , \hspace3mm \Big|\frac{i-1}{n} - U_{(i)} \Big| \Big) \] with $U_i \sim_{iid} Unif(0,1)$.
No matter what $F$ is, distribution of the test statistic KS under the null is known.
KS test with completely specified $F_0$ is a distribution-free test.
When $F_0$ contains nuisance parameter, Distribution of $D_n$ under the null depends on $F$. KS test becomes only asymptotically distribution free test. For finite $n$, null distribution must be computed by Monte Carlo simulation.

K-S test in Time Series Setting

In time sries, we want to use K-S test for checking distribution of innovations. (errors).
Since we assume $e_t$ are i.i.d. noise, if we can observe $e_t$, K-S test can be directly applicable.
However, in Time Series analysis $e_t$ are not observable, and only residuals $\hat e_t$ are available.
Can we still use K-S test on $\hat e_t$?

  Fn <- ecdf(as.numeric(res1))   #- Fn becomes an function
  plot( ecdf(as.numeric(res1)) )

  t=seq(-4,4, .01)
  plot(t,  Fn(t), type="l" )
  lines(t, pnorm(t,0,1), col="red")

Lecture 27 - ARMA-GARCH model

Nao Mimoto - Dept. of Statistics : The University of Akron

Apr 4 2018